{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This notebook shows how to make a 3D geological model that you can then be embedded in a website"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a notebook showing how to recreate a model of a geothermal reservoir using publicly available data.\n",
    "The model is of a reservoir in Utah that is part of a project called FORGE (Frontier Observatory for Research in Geothermal Energy).\n",
    "\n",
    "See image that you will generate here: https://github.com/ahinoamp/Example3DGeologicModelUsingVTKPlotter/blob/master/ModelImage.png\n",
    "\n",
    "Original data source links are shown in the end."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Image of Yaktocat](https://raw.githubusercontent.com/ahinoamp/Example3DGeologicModelUsingVTKPlotter/master/ModelImage.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#imports\n",
    "from vedo import *\n",
    "import pandas as pd\n",
    "from scipy.spatial import Delaunay\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import numpy as np\n",
    "    \n",
    "embedWindow('k3d') # panel or itkwidgets or False (for a popup)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import the file from github\n",
    "basegithuburl = 'https://raw.githubusercontent.com/ahinoamp/Example3DGeologicModelUsingVTKPlotter/master/'\n",
    "    \n",
    "#Load surfaces\n",
    "fileVertices = basegithuburl+'land_surface_vertices.csv'\n",
    "landSurfacePD =pd.read_csv(fileVertices)\n",
    "\n",
    "fileVertices = basegithuburl+'175C_vertices.csv'\n",
    "vertices_175CPD =pd.read_csv(fileVertices)\n",
    "\n",
    "fileVertices = basegithuburl+'225C_vertices.csv'\n",
    "vertices_225CPD =pd.read_csv(fileVertices)\n",
    "\n",
    "fileVertices = basegithuburl+'Negro_Mag_Fault_vertices.csv'\n",
    "Negro_Mag_Fault_verticesPD =pd.read_csv(fileVertices)\n",
    "\n",
    "fileVertices = basegithuburl+'Opal_Mound_Fault_vertices.csv'\n",
    "Opal_Mound_Fault_verticesPD =pd.read_csv(fileVertices)\n",
    "\n",
    "fileVertices = basegithuburl+'top_granitoid_vertices.csv'\n",
    "top_granitoid_verticesPD =pd.read_csv(fileVertices)\n",
    "\n",
    "fileVertices = basegithuburl+'top_granitoid_vertices.csv'\n",
    "border =pd.read_csv(fileVertices)\n",
    "\n",
    "fileVertices = basegithuburl+'Microseismic.csv'\n",
    "microseismic =pd.read_csv(fileVertices)\n",
    "\n",
    "#The well path and different logs for the well paths\n",
    "filepath = basegithuburl+'path5832.csv'\n",
    "well_5832_path =pd.read_csv(filepath)\n",
    "\n",
    "filepath = basegithuburl+'temperature5832.csv'\n",
    "temp_well =pd.read_csv(filepath)\n",
    "\n",
    "filepath = basegithuburl+'nphi5832.csv'\n",
    "nphi_well =pd.read_csv(filepath)\n",
    "\n",
    "filepath = basegithuburl+'pressure5832.csv'\n",
    "pressure_well =pd.read_csv(filepath)\n",
    "\n",
    "#Since most of the wells in the area were just vertical, I split them into two files:\n",
    "#One file with the top of the wells and the other with the bottom point of the wellbore\n",
    "file = basegithuburl + 'MinPointsWells.csv'\n",
    "wellsmin =pd.read_csv(file)\n",
    "file = basegithuburl + 'MaxPointsWells.csv'\n",
    "wellsmax =pd.read_csv(file)\n",
    "\n",
    "#Project boundary area on the surface\n",
    "file = basegithuburl + 'FORGE_Border.csv'\n",
    "border = pd.read_csv(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a plot\n",
    "plot = Plotter(axes=1, bg='white', interactive=1)\n",
    "\n",
    "####################\n",
    "## 1. land surface: a mesh with varying color\n",
    "####################\n",
    "\n",
    "#perform a 2D Delaunay triangulation to get the cells from the point cloud\n",
    "tri = Delaunay(landSurfacePD.values[:, 0:2])\n",
    "#create a mesh object for the land surface\n",
    "landSurface = Mesh([landSurfacePD.values, tri.simplices])\n",
    "\n",
    "#in order to color it by the elevation, we extract the z value\n",
    "elevation = landSurface.cellCenters()[:, 2]   # pick z coordinates of cells\n",
    "\n",
    "#unfortunatly I couldn't find a good colormap for terrain without ocean. \n",
    "#So we'll need to truncate the \"terrain\" color map\n",
    "cmap = plt.get_cmap('terrain')\n",
    "truncateMin = 0.23 #We want to start the colormap about a quarter of the way in until the end\n",
    "truncateMax = 1\n",
    "cmap_terrain_no_ocean = colors.LinearSegmentedColormap.from_list(\n",
    "    'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=truncateMin, b=truncateMax),\n",
    "    cmap(np.linspace(truncateMin, truncateMax, 100)))\n",
    "\n",
    "#Now we color the cell by elevation\n",
    "landSurface.cmap(cmap_terrain_no_ocean, elevation, on='cells')\n",
    "\n",
    "#We give the object a name\n",
    "landSurface.name='Land Surface'\n",
    "\n",
    "#We add it to the plot\n",
    "plot += landSurface\n",
    "\n",
    "\n",
    "####################\n",
    "## 2. Different meshes with constant colors\n",
    "####################\n",
    "\n",
    "#Mesh of 175 C isotherm\n",
    "tri = Delaunay(vertices_175CPD.values[:, 0:2])\n",
    "vertices_175C = Mesh([vertices_175CPD.values, tri.simplices]).c(\"orange\").opacity(0.3)\n",
    "vertices_175C.name='175C temperature isosurface'\n",
    "plot += vertices_175C\n",
    "\n",
    "#Mesh of 225 C isotherm\n",
    "tri = Delaunay(vertices_225CPD.values[:, 0:2])\n",
    "vertices_225CT = Mesh([vertices_225CPD.values, tri.simplices]).c(\"red\").opacity(0.4)\n",
    "vertices_225CT.name='225C temperature isosurface'\n",
    "plot += vertices_225CT\n",
    "\n",
    "#Negro fault\n",
    "tri = Delaunay(Negro_Mag_Fault_verticesPD.values[:, 1:3])\n",
    "Negro_Mag_Fault_vertices = Mesh([Negro_Mag_Fault_verticesPD.values, tri.simplices]).c(\"f\").opacity(0.4)\n",
    "Negro_Mag_Fault_vertices.name='Negro Fault'\n",
    "plot += Negro_Mag_Fault_vertices\n",
    "\n",
    "#Opal fault\n",
    "tri = Delaunay(Opal_Mound_Fault_verticesPD.values[:, 1:3])\n",
    "Opal_Mound_Fault_vertices = Mesh([Opal_Mound_Fault_verticesPD.values, tri.simplices]).c(\"g\").opacity(0.4)\n",
    "Opal_Mound_Fault_vertices.name='Opal Mound Fault'\n",
    "plot += Opal_Mound_Fault_vertices\n",
    "\n",
    "#Top Granite\n",
    "xyz = top_granitoid_verticesPD.values\n",
    "xyz[:, 2] = top_granitoid_verticesPD.values[:,2]-20\n",
    "tri = Delaunay(top_granitoid_verticesPD.values[:, 0:2])\n",
    "top_granitoid_vertices = Mesh([xyz, tri.simplices]).c(\"darkcyan\")\n",
    "top_granitoid_vertices.name='Top of granite surface'\n",
    "plot += top_granitoid_vertices\n",
    "\n",
    "####################\n",
    "## 3. Point objects\n",
    "####################\n",
    "\n",
    "#FORGE Boundary\n",
    "#Since the boundary area did not have a Z column, I assigned a Z value for where I wanted it to appear\n",
    "border['zcoord'] = 1650\n",
    "borderxyz = border[['xcoord', 'ycoord', 'zcoord']]\n",
    "boundary = Points(borderxyz.values).c('k')\n",
    "boundary.name='FORGE area boundary'\n",
    "plot+=boundary\n",
    "\n",
    "#Microseismic\n",
    "microseismicxyz = microseismic[['xloc','yloc','zloc']]\n",
    "scals = microseismic[['mw']]\n",
    "microseismicPts = Points(microseismicxyz.values, r=3).cmap(\"jet\", scals)\n",
    "microseismicPts.name='Microseismic events'\n",
    "plot+=microseismicPts\n",
    "\n",
    "####################\n",
    "## 4. Line objects\n",
    "####################\n",
    "\n",
    "#The path of well 58_32\n",
    "xyz = well_5832_path[['X', 'Y', 'Z']].values\n",
    "Well = Line(xyz)\n",
    "Well.name='Well 58-32'\n",
    "plot+=Well\n",
    "\n",
    "#A porosity log in the well\n",
    "xyz = nphi_well[['X', 'Y', 'Z']].values\n",
    "porosity = nphi_well['Nphi'].values\n",
    "Well = Line(xyz).c('gold')\n",
    "Well.name='Porosity log well 58-32'\n",
    "plot+=Well\n",
    "\n",
    "#This well data is actually represented by points since as of right now, \n",
    "#since the k3d embedding does not support colors on the lines, and I wanted to show the colors\n",
    "xyz = pressure_well[['X', 'Y', 'Z']].values\n",
    "pressure = pressure_well['Pressure'].values\n",
    "Well = Points(xyz, r=1).cmap(\"cool\", pressure)\n",
    "Well.name='Pressure log well 58-32'\n",
    "plot+=Well\n",
    "\n",
    "#Temperatue log\n",
    "xyz = temp_well[['X', 'Y', 'Z']].values\n",
    "scals = temp_well['Temperature'].values\n",
    "Well = Points(xyz, r=1).cmap(\"seismic\", scals)\n",
    "Well.name='Temperature log well 58-32'\n",
    "plot+=Well\n",
    "\n",
    "\n",
    "####################\n",
    "## 5. Multi-line objects\n",
    "####################\n",
    "\n",
    "#There is some preprocessing that needs to be done here in order to get two lists of points\n",
    "#defining the start and end of the lines that will be representing the wellbores\n",
    "xyzmin = wellsmin[['x', 'y', 'z']].values\n",
    "xyzmax = wellsmax[['x', 'y', 'z']].values\n",
    "p1=[]\n",
    "p2=[]\n",
    "for i in range(len(xyzmin)):\n",
    "    p1.append(xyzmin[i,:])\n",
    "    p2.append(xyzmax[i,:])\n",
    "\n",
    "Wells = Lines(p1, p2, c='gray', alpha=1, lw=3)\n",
    "Wells.name='Pre-existing wellbores'\n",
    "plot+=Wells\n",
    "\n",
    "####################\n",
    "## 6. Done. show the plot\n",
    "####################\n",
    "plot.show(viewup='z')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To embed this model in a website, you have to be in \"k3d\" embedding mode, click on \"controls\", and then on \"Snapshot HTML\". This should prompt you to save your HTML file on your system! You can now embed that HTML code in a container on your website! yay!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the data is from the Geothermal Data Repository (https://gdr.openei.org/home) uploaded by the Energy and Geoscience Institute at the University of Utah.\n",
    "\n",
    "1. Earthquake data: \"Utah FORGE: Earthquake Catalog\", https://gdr.openei.org/submissions/1039\n",
    "2. Well 58-32  porosity log: \"Utah FORGE: Well 58-32 Schlumberger FMI Logs DLIS and XML files\", https://gdr.openei.org/submissions/1076\n",
    "3. Well 58-32 pressure and temperature logs: \"Utah FORGE: Milford Deep Test Well 58-32 (MU-ESW1) Pressure and Temperature Logs\", https://gdr.openei.org/submissions/1101\n",
    "4. Microseismic data: \"Utah FORGE: Microseismic Events\", https://gdr.openei.org/submissions/1151\n",
    "5. Well data and surfaces data: \"Utah FORGE: Well Data for Student Competition\", https://gdr.openei.org/submissions/1111"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
